library(meta)
library(metafor)
library(ggrepel)
library(dmetar)
library(smd)
library(tidyverse)
library(PerformanceAnalytics)
library(ggplot2); theme_set(theme_bw(base_size = 15)+
theme(axis.text.x = element_text(angle = 45, hjust=1),
strip.text = element_text(colour = 'black', face="bold",size=12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(size = 0.7),
axis.ticks.length=unit(.10, "cm"),
axis.ticks = element_line(size=0.7),
strip.background = element_blank()))
##Loading data and prepping tables It is worth noting that using MD results in higher heterogeneity values than using SMD. Here I use MD as we have discussed. The significance does not change, but it affects the multi-model inference later. Using SMD, there is no significant influence… using MD there is.
Here we see that MOG papers are biased. We will have to state this as a limitation
# Generate funnel plot (we do not include study labels here)
funnel.meta(mogmetagen)
title("MOG ARR")
# Plot over points
eggers.test(mogmetagen)
## Eggers' test of the intercept
## =============================
##
## intercept 95% CI t p
## -2.654 -3.62 - -1.69 -5.404 5.847195e-05
##
## Eggers' test indicates the presence of funnel plot asymmetry.
NMOSD papers are much more even.
# Generate funnel plot (we do not include study labels here)
funnel.meta(nmosdmetagen)
title("NMOSD ARR")
# Plot over points
eggers.test(nmosdmetagen)
## Warning in metabias.meta(x, k.min = 3, method = "linreg"): 1 observation(s)
## dropped due to missing values
## Eggers' test of the intercept
## =============================
##
## intercept 95% CI t p
## -0.591 -2.2 - 1.02 -0.72 0.4797325
##
## Eggers' test does not indicate the presence of funnel plot asymmetry.
MOG EDSS scores are not biased across studies.
# Generate funnel plot (we do not include study labels here)
funnel.meta(mogedssmetagen)
title("MOG EDSS")
# Plot over points
eggers.test(mogedssmetagen)
## Warning in metabias.meta(x, k.min = 3, method = "linreg"): 1 observation(s)
## dropped due to missing values
## Eggers' test of the intercept
## =============================
##
## intercept 95% CI t p
## -0.265 -2 - 1.47 -0.298 0.7730529
##
## Eggers' test does not indicate the presence of funnel plot asymmetry.
NMOSD EDSS scores are also not biased.
# Generate funnel plot (we do not include study labels here)
funnel.meta(nmosdedssmetagen)
title("NMOSD EDSS")
# Plot over points
eggers.test(nmosdedssmetagen)
## Warning in metabias.meta(x, k.min = 3, method = "linreg"): 1 observation(s)
## dropped due to missing values
## Eggers' test of the intercept
## =============================
##
## intercept 95% CI t p
## -0.257 -1.48 - 0.97 -0.411 0.6878221
##
## Eggers' test does not indicate the presence of funnel plot asymmetry.
Here a random effects model reveals significant difference in mean ARR with lower ARRs post RTX.
##plot mog ARR forest plot
forest.meta(mogmetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
leftlabs = c("Study", "g", "SE"))
print(mogmetagen)
## Review: Mogmeta
##
## Number of studies combined: k = 18
##
## MD 95%-CI z p-value
## Random effects model -0.9213 [-1.2424; -0.6002] -5.62 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.2650 [0.0813; 0.9016]; tau = 0.5148 [0.2851; 0.9495]
## I^2 = 87.2% [81.3%; 91.3%]; H = 2.80 [2.31; 3.39]
##
## Test of heterogeneity:
## Q d.f. p-value
## 133.01 17 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
As shown before in Damata et al., NMOSD post RTX ARRs are significantly lower.
##plot mog ARR forest plot
forest.meta(nmosdmetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
leftlabs = c("Study", "g", "SE"))
print(nmosdmetagen)
## Review: NMOSDmeta
##
## Number of studies combined: k = 23
##
## MD 95%-CI z p-value
## Random effects model -1.7324 [-2.1507; -1.3141] -8.12 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.7079 [0.3618; 2.6610]; tau = 0.8414 [0.6015; 1.6313]
## I^2 = 79.0% [69.1%; 85.7%]; H = 2.18 [1.80; 2.65]
##
## Test of heterogeneity:
## Q d.f. p-value
## 104.64 22 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
The merge analysis reveals that NMOSD patients benefit more than MOG patients. Between groups Q=9.09 p=0.003.
##plot mog and nmosd together
mergemetagen <- update.meta(mergemetagen,
subgroup = Group,
tau.common = FALSE)
forest.meta(mergemetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
subgroup = TRUE,
leftlabs = c("Study", "MD", "SE"))
print(mergemetagen)
## Review: mergemeta
##
## Number of studies combined: k = 41
##
## MD 95%-CI z p-value
## Random effects model -1.3866 [-1.6791; -1.0942] -9.29 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.6122 [0.3644; 1.4990]; tau = 0.7824 [0.6037; 1.2243]
## I^2 = 93.2% [91.7%; 94.5%]; H = 3.85 [3.47; 4.27]
##
## Test of heterogeneity:
## Q d.f. p-value
## 592.56 40 < 0.0001
##
## Results for subgroups (random effects model):
## k MD 95%-CI tau^2 tau Q I^2
## Group = MOG 18 -0.9213 [-1.2424; -0.6002] 0.2650 0.5148 133.01 87.2%
## Group = NMOSD 23 -1.7324 [-2.1507; -1.3141] 0.7079 0.8414 104.64 79.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 9.09 1 0.0026
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
EDSS appeared to be significantly better in those post RTX.
##plot MOG EDSS
forest.meta(mogedssmetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
leftlabs = c("Study", "g", "SE"))
print(mogedssmetagen)
## Review: MogEDSSmeta
##
## Number of studies combined: k = 10
##
## MD 95%-CI z p-value
## Random effects model -0.8384 [-1.4146; -0.2621] -2.85 0.0044
##
## Quantifying heterogeneity:
## tau^2 = 0.4549 [0.0231; 3.3570]; tau = 0.6744 [0.1520; 1.8322]
## I^2 = 55.4% [9.4%; 78.1%]; H = 1.50 [1.05; 2.14]
##
## Test of heterogeneity:
## Q d.f. p-value
## 20.20 9 0.0167
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
EDSS is significantly improved post RTX in NMOSD patients.
##plot MOG EDSS
forest.meta(nmosdedssmetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
leftlabs = c("Study", "g", "SE"))
print(nmosdedssmetagen)
## Review: NMOSDEDSSmeta
##
## Number of studies combined: k = 15
##
## MD 95%-CI z p-value
## Random effects model -1.1439 [-1.6257; -0.6620] -4.65 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.2571 [0.0000; 1.5773]; tau = 0.5071 [0.0000; 1.2559]
## I^2 = 29.4% [0.0%; 62.0%]; H = 1.19 [1.00; 1.62]
##
## Test of heterogeneity:
## Q d.f. p-value
## 19.83 14 0.1357
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
Here we see that there is no difference between EDSS scores of patients pre and post RTX between MOG and NMOSD (p=0.4254)
##plot mog and nmosd together
mergeedssmetagen <- update.meta(mergeedssmetagen,
subgroup = Group,
tau.common = FALSE)
forest.meta(mergeedssmetagen,
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
subgroup = TRUE,
leftlabs = c("Study", "MD", "SE"))
print(mergeedssmetagen)
## Review: mergeEDSSmeta
##
## Number of studies combined: k = 25
##
## MD 95%-CI z p-value
## Random effects model -0.9903 [-1.3590; -0.6217] -5.27 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.3443 [0.0215; 1.2263]; tau = 0.5868 [0.1468; 1.1074]
## I^2 = 42.2% [6.9%; 64.1%]; H = 1.32 [1.04; 1.67]
##
## Test of heterogeneity:
## Q d.f. p-value
## 41.51 24 0.0146
##
## Results for subgroups (random effects model):
## k MD 95%-CI tau^2 tau Q I^2
## Group = MOG 10 -0.8384 [-1.4146; -0.2621] 0.4549 0.6744 20.20 55.4%
## Group = NMOSD 15 -1.1439 [-1.6257; -0.6620] 0.2571 0.5071 19.83 29.4%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 0.64 1 0.4254
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
Here we removed mean baseline EDSS and prop children. As they significantly correlate with other predictors.
##meta-regression
mogreg<- read.csv("/Users/bosun/Documents/DPhil/OANG/6\ -\ Collaborations/2021-12-27_Val_meta/data/2022-03-27-mog_meta-regression data_last.csv")
mogreg <- mogreg %>%
dplyr::left_join(mogmeta, by = c("STUDY", "n"))
mogreg$weights <- 1/sqrt(mogreg$vi)
plot_metareg <- function(data = mogreg, variable = NULL){
mogreg2 <- mogreg[!is.na(mogreg[[paste0(variable)]]),]
reg <- rma(yi = yi,
vi = vi,
mods = mogreg2[[paste0(variable)]],
data = mogreg2)
# Specify basic plot, mapping sex to the x-axis, effect size 'd' to the y-axis,
# and 'weights' to the weight parameter.
preds <- predict.rma(reg)
plotdb <- cbind(mogreg2, preds)
ggplot(plotdb, aes(x = plotdb[[paste0(variable)]], y = plotdb[["yi"]], label=plotdb[["STUDY"]])) +
geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha=0.05)+
geom_line(data = plotdb,aes(x = plotdb[[paste0(variable)]],y = plotdb$pred, size = 0.5), linetype="dashed", alpha=0.1)+
geom_text_repel(min.segment.length = 0, seed = 42, box.padding = 0.5, size=3, force=1)+
geom_point(aes(size=weights),shape = 21, alpha=1, fill="white") +
theme(legend.position = "none")+
xlab(variable)+
ylab("MD")
}
#Get modifiers
plotparam<- c("mean_time_to_RTX_months",
"mean_onset_age",
"mean_fup_months",
"mean_baseline_ARR",
"mean_baseline_EDSS",
"prop_relapsing",
"prop_children",
"prop_female")
#filter modifiers for highly correlated
mogreg[,unlist(plotparam)] %>%
chart.Correlation()
##remove prop_children as correlated with mean age onset
## remove mean EDSS as correlated highly with age of onset
plotparam <-plotparam[!plotparam %in% c("mean_baseline_EDSS", "prop_children")]
reg_plotlist <- lapply(c(plotparam,"n"), function(x){
p <- plot_metareg(data=mogreg, variable = x)
return(p)
})
##plot bubble plots
patchwork::wrap_plots(reg_plotlist)
##get univariate model metrics
get_unimodels <- function(data, mod){
m_reg <- rma.uni(yi = yi,
vi = vi,
mods = data[[paste0(mod)]],
data = data)
metrics <- c(mod = mod, tau2 = round(m_reg$tau2,3),
tau =round(sqrt(m_reg$tau2), 3),
I2 = round(m_reg$I2, 3),
R2 = round(m_reg$R2, 3),
beta = round(m_reg$beta["mods",],3),
zval = round(m_reg$zval[2],3),
pval = round(m_reg$pval[2],7),
ci.lb = round(m_reg$ci.lb[2],3),
ci.ub = round(m_reg$ci.ub[2],3))
}
modregs <- lapply(c(plotparam,"n"), function(x){
get_unimodels(mogreg, x)
})
modregs <- arrange(as.data.frame(do.call(rbind, modregs)), desc(as.numeric(R2)))
kableExtra::kable(modregs)%>% kableExtra::kable_material(c("striped", "hover"))
| mod | tau2 | tau | I2 | R2 | beta.mods | zval | pval | ci.lb | ci.ub |
|---|---|---|---|---|---|---|---|---|---|
| mean_time_to_RTX_months | 0 | 0 | 0 | 100 | 0.016 | 4.118 | 3.82e-05 | 0.008 | 0.024 |
| mean_baseline_ARR | 0 | 0 | 0 | 100 | -0.85 | -5.035 | 5e-07 | -1.18 | -0.519 |
| prop_female | 0.196 | 0.443 | 53.192 | 50.714 | -3.67 | -1.922 | 0.0546228 | -7.413 | 0.073 |
| n | 0.194 | 0.441 | 48.954 | 10.659 | -0.006 | -0.954 | 0.3402878 | -0.017 | 0.006 |
| mean_onset_age | 0.303 | 0.55 | 59.152 | 0 | 0.028 | 1.519 | 0.1286385 | -0.008 | 0.063 |
| mean_fup_months | 0.239 | 0.489 | 59.02 | 0 | -0.008 | -0.651 | 0.5151981 | -0.03 | 0.015 |
| prop_relapsing | 0.244 | 0.494 | 58.782 | 0 | -0.978 | -0.692 | 0.4889695 | -3.747 | 1.791 |
Here, all possible predictor combinations with multi-model inference, we see that mean ARR is the most important predictor of RTX response.
multimodel.inference(TE = "yi",
seTE = "sei",
data = mogreg,
predictors = c(plotparam),
interaction = TRUE)
## You entered 6 predictors. Interactions can only be modeled for four or less predictors. Therefore, no interactions are modeled.
##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 6%
|
|===== | 8%
|
|======= | 9%
|
|======== | 11%
|
|========= | 12%
|
|========== | 14%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 23%
|
|================== | 25%
|
|=================== | 27%
|
|==================== | 28%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 41%
|
|============================== | 42%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 58%
|
|========================================== | 59%
|
|=========================================== | 61%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 86%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 91%
|
|================================================================= | 92%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
##
##
## Multimodel Inference: Final Results
## --------------------------
##
## - Number of fitted models: 64
## - Full formula: ~ mean_time_to_RTX_months + mean_onset_age + mean_fup_months + mean_baseline_ARR + prop_relapsing + prop_female
## - Coefficient significance test: knha
## - Interactions modeled: no
## - Evaluation criterion: AICc
##
##
## Best 5 Models
## --------------------------
##
##
## Global model call: metafor::rma(yi = TE, sei = seTE, mods = form, data = glm.data,
## method = method, test = test)
## ---
## Model selection table
## (Int) men_bsl_ARR men_fup_mnt men_ons_age prp_fml prp_rlp df logLik AICc
## 2 + -0.8497 3 -5.647 19.7
## 18 + -0.7958 -1.907 4 -2.903 21.8
## 4 + -0.8756 0.006311 4 -4.954 22.4
## 6 + -0.8607 0.01268 4 -4.452 22.6
## 34 + -0.8387 -0.7559 4 -5.177 22.8
## delta weight
## 2 0.00 0.486
## 18 2.11 0.169
## 4 2.66 0.129
## 6 2.92 0.113
## 34 3.11 0.103
## Models ranked by AICc(x)
##
##
## Multimodel Inference Coefficients
## --------------------------
##
##
## Estimate Std. Error z value Pr(>|z|)
## intrcpt 0.3983073530 0.602707328 0.66086363 0.5086998
## mean_baseline_ARR -0.8381088054 0.150244981 5.57828156 0.0000000
## prop_female -0.3178070621 0.800662121 0.39693031 0.6914189
## mean_fup_months 0.0008785040 0.002877126 0.30534082 0.7601066
## mean_onset_age 0.0016179293 0.005125106 0.31568701 0.7522401
## prop_relapsing -0.0923513767 0.348598158 0.26492216 0.7910694
## mean_time_to_RTX_months 0.0001573002 0.001721205 0.09138954 0.9271831
##
##
## Predictor Importance
## --------------------------
##
##
## model importance
## 1 mean_baseline_ARR 0.98460599
## 2 prop_female 0.16270696
## 3 mean_fup_months 0.14308581
## 4 mean_onset_age 0.12012746
## 5 prop_relapsing 0.11775970
## 6 mean_time_to_RTX_months 0.08825281
## $relapsing_mog
## Review: relapsing_mog
##
## Number of studies combined: k = 16
##
## MD 95%-CI z p-value
## Random effects model -1.0270 [-1.3226; -0.7314] -6.81 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1369 [0.0045; 0.8841]; tau = 0.3700 [0.0673; 0.9403]
## I^2 = 47.6% [6.3%; 70.7%]; H = 1.38 [1.03; 1.85]
##
## Test of heterogeneity:
## Q d.f. p-value
## 28.60 15 0.0181
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
##
## $adult_mog
## Review: adult_mog
##
## Number of studies combined: k = 10
##
## MD 95%-CI z p-value
## Random effects model -1.0018 [-1.5620; -0.4416] -3.50 0.0005
##
## Quantifying heterogeneity:
## tau^2 = 0.4621 [0.0956; 1.4413]; tau = 0.6798 [0.3092; 1.2006]
## I^2 = 85.5% [75.2%; 91.6%]; H = 2.63 [2.01; 3.44]
##
## Test of heterogeneity:
## Q d.f. p-value
## 62.27 9 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
##
## $children_mog
## Review: children_mog
##
## Number of studies combined: k = 7
##
## MD 95%-CI z p-value
## Random effects model -1.1957 [-1.8592; -0.5322] -3.53 0.0004
##
## Quantifying heterogeneity:
## tau^2 = 0.5047 [0.0592; 6.4154]; tau = 0.7105 [0.2432; 2.5329]
## I^2 = 65.2% [21.7%; 84.5%]; H = 1.69 [1.13; 2.54]
##
## Test of heterogeneity:
## Q d.f. p-value
## 17.23 6 0.0085
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
##
## $fu12_mog
## Review: fu12_mog
##
## Number of studies combined: k = 14
##
## MD 95%-CI z p-value
## Random effects model -1.1253 [-1.5243; -0.7264] -5.53 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.2533 [0.0100; 1.2335]; tau = 0.5033 [0.0998; 1.1106]
## I^2 = 50.3% [8.2%; 73.1%]; H = 1.42 [1.04; 1.93]
##
## Test of heterogeneity:
## Q d.f. p-value
## 26.15 13 0.0162
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
##
## $first_mog
## Review: firstline_mog
##
## Number of studies combined: k = 5
##
## MD 95%-CI z p-value
## Random effects model -1.3202 [-2.6540; 0.0136] -1.94 0.0524
##
## Quantifying heterogeneity:
## tau^2 = 1.8373 [0.4004; 17.6857]; tau = 1.3555 [0.6327; 4.2054]
## I^2 = 84.9% [66.3%; 93.2%]; H = 2.57 [1.72; 3.84]
##
## Test of heterogeneity:
## Q d.f. p-value
## 26.44 4 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
for (i in 1:length(subgroupsmetagen)){
forest.meta(subgroupsmetagen[[paste0(names(subgroupsmetagen[i]))]],
sortvar = yi,
predict = TRUE,
print.tau2 = FALSE,
leftlabs = c("Study", "g", "SE"))
grid::grid.text(paste0(names(subgroupsmetagen[i])), .5, .9, gp=grid::gpar(cex=2))
}